{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import fst\n",
    "import matplotlib.pyplot as P\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# HMM简介\n",
    "\n",
    "隐马尔可夫模型(HMMs)是一个用于建模时序数据的模型。它的优点是非常直观而且容易实现。你可以参考[Wikipedia](https://en.wikipedia.org/wiki/Hidden_Markov_model) 或者google搜索到很多相关资料。([这篇文章](http://www.comp.leeds.ac.uk/roger/HiddenMarkovModels/html_dev/main.html) 是我最喜欢的文章之一。\n",
    "\n",
    "对于语音识别来说，首先需要阅读的是[Rabiner]的论文(http://www.ee.columbia.edu/~dpwe/e6820/papers/Rabiner89-hmm.pdf)，这篇论文非常写到很容易理解而且是专门讨论它在语音识别中的应用。如果这是你第一次了解它，你可能需要先阅读这些相关资料。这个notebook会简单的回顾它的一些要点然后展示怎么用Python来实现。\n",
    "\n",
    "\n",
    "## HMM 算法\n",
    "\n",
    "Rabiner的论文提到了HMM有3个问题： \n",
    "\n",
    "  1. 给定模型(参数)，计算观察的概率，也就是P(O|M) - 这需要使用前向算法 \n",
    "  2. 给定模型(参数)和观察，计算最可能的状态序列，也就是$\\underset{S}{argmax}P(S|O;M)$。这需要使用Viterbi算法。 \n",
    "  3. 学习问题，调整模型参数，使得概率P(O|M)最大(最大似然)。这需要Baum-Welch算法，也叫前向后向算法。\n",
    "  \n",
    "HMM有如下的参数定义：\n",
    "\n",
    "  * 发射概率 - 通常我们使用高斯混合模型(GMM)，但是我们也可以使用深度神经网络(比如MLP)。为了与Rabiner的论文一致，我们假设输出是离散的。因此可以用概率分布函数$b_{jk}$来描述。它表示给定状态j出现观察k的概率。\n",
    "  * 隐状态 - 在我们的例子中，隐状态是我们需要识别的单元：比如因素(phonemes)，但是也有可能是其它的单元比如词。 \n",
    "  * 转移概率 - 从一个状态跳转到另一个状态的概率  \n",
    "  * 先验概率 - 系统处于不同初始状态的概率。(如果忽略也问题不大)\n",
    "  \n",
    "为了便于阅读，我们这里使用和Rabiner论文一致的记号:\n",
    "  * 发射概率我们用$b_{jk}$表示。\n",
    "  * 跳转概率用$a_{ij}$表示。\n",
    "  * 先验概率用$pi_i$表示。\n",
    "  \n",
    "## 示例模型\n",
    "\n",
    "这里我们会随机的生成一些值来演示算法。我们这里的示例HMM有3个状态，观察有4种可能值。观察的长度是5。\n",
    "\n",
    "我们把它们定义成常量。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "N=3\n",
    "M=4\n",
    "T=5"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 观察序列\n",
    "\n",
    "我们先随机生成一个观察序列："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[2 2 0 2 1]\n"
     ]
    }
   ],
   "source": [
    "O=np.random.randint(0,M,T)\n",
    "print(O)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "发射概率可以用一个矩阵来表示。这个矩阵的大小是3x4，每一行表示一个状态，$b_{jk}$表示状态j发射观察k的概率。\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[0.48229459 0.27721979 0.12239203 0.11809359]\n",
      " [0.11822906 0.14038183 0.42902164 0.31236747]\n",
      " [0.34614273 0.10475221 0.38320261 0.16590244]]\n"
     ]
    }
   ],
   "source": [
    "b=np.random.random((N,M))\n",
    "b/=np.array(b.sum(axis=1, keepdims=True))\n",
    "print(b)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "因为是概率，因此要求$\\sum_kb_{jk}=1$，也就是每一行加起来等于1。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1. 1. 1.]\n"
     ]
    }
   ],
   "source": [
    "print(b.sum(axis=1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "给定这些概率和一个观察序列，我们很容易计算的不同时刻处于不同状态的发射概率："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[0.12239203 0.42902164 0.38320261]\n",
      " [0.12239203 0.42902164 0.38320261]\n",
      " [0.48229459 0.11822906 0.34614273]\n",
      " [0.12239203 0.42902164 0.38320261]\n",
      " [0.27721979 0.14038183 0.10475221]]\n"
     ]
    }
   ],
   "source": [
    "ob=np.zeros((T,N))\n",
    "for t in range(T):\n",
    "    for i in range(N):\n",
    "        ob[t,i]=b[i,O[t]]\n",
    "print(ob)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "或者更加简洁高效的写法："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[0.12239203 0.42902164 0.38320261]\n",
      " [0.12239203 0.42902164 0.38320261]\n",
      " [0.48229459 0.11822906 0.34614273]\n",
      " [0.12239203 0.42902164 0.38320261]\n",
      " [0.27721979 0.14038183 0.10475221]]\n"
     ]
    }
   ],
   "source": [
    "ob=b.T[O]\n",
    "print(ob)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "注意我们这里的示例使用的是离散的输出，我们用一个发射概率矩阵就可以表示了，但是在实际的系统中，我们通常使用一个模型，比如GMM或者MLP来表示这个发射概率。这里为了简单起见我们使用了离散的输出，但是读者应该知道在实际的系统中很少使用离散的输出。\n",
    "\n",
    "#### 状态转移矩阵\n",
    "\n",
    "状态转移矩阵是一个3x3的矩阵，它的每一行加起来是1："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[0.39884173 0.3486846  0.25247367]\n",
      " [0.08072832 0.59021607 0.32905562]\n",
      " [0.25993654 0.25873506 0.4813284 ]]\n",
      "[1. 1. 1.]\n"
     ]
    }
   ],
   "source": [
    "a=np.random.random((N,N))\n",
    "a/=np.array([a.sum(axis=-1)]).T\n",
    "print(a)\n",
    "print(a.sum(axis=1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 先验概率\n",
    "\n",
    "下面我们随机的生成先验概率。我们可以设置随机值，也可以让所有的状态的先验概率相等。\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0.23738012 0.47646268 0.28615719]\n"
     ]
    }
   ],
   "source": [
    "pi=np.random.random(N)\n",
    "pi/=pi.sum()\n",
    "print(pi)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 最终的模型\n",
    "\n",
    "为了更容易阅读，我们把这些概率放到一个tuple里："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "hmm=(O,pi,a,b,N,M,T)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "下面把模型保存到磁盘："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pickle\n",
    "\n",
    "with open('../data/hmm.pkl','wb') as f:\n",
    "    pickle.dump(hmm,f,pickle.HIGHEST_PROTOCOL)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 前向算法\n",
    "\n",
    "下面会直接实现Rabiner的论文里的算法。注意原文的下标是1开始的，而程序的下标是0开始。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[1.05337157e-01 1.57878382e-01 6.43555966e-02]\n",
      " [4.00636322e-02 4.24970349e-02 2.45214076e-02]\n",
      " [8.41601229e-03 2.33586277e-03 2.20264128e-02]\n",
      " [1.66340860e-03 5.68899442e-03 2.66746419e-03]\n",
      " [9.48980897e-04 1.19242838e-04 2.32921745e-03]]\n"
     ]
    }
   ],
   "source": [
    "def forward(O,pi,a,b,N,M,T):\n",
    "    fwd=np.zeros((T,N))\n",
    "\n",
    "    #初始化\n",
    "    for i in range(N):\n",
    "        fwd[0,i]=pi[i]*b[i,O[0]]\n",
    "\n",
    "    #递推\n",
    "    for t in range(T-1):\n",
    "        for j in range(N):\n",
    "            s=0\n",
    "            for i in range(N):\n",
    "                s+=fwd[t,i]*a[i,j]\n",
    "            fwd[t+1,j]=s*b[j,O[t+1]]\n",
    "            \n",
    "    return fwd\n",
    "\n",
    "print(forward(*hmm))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "我们也可以使用Numpy更加高效的实现："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[1.05337157e-01 1.57878382e-01 6.43555966e-02]\n",
      " [4.00636322e-02 4.24970349e-02 2.45214076e-02]\n",
      " [8.41601229e-03 2.33586277e-03 2.20264128e-02]\n",
      " [1.66340860e-03 5.68899442e-03 2.66746419e-03]\n",
      " [9.48980897e-04 1.19242838e-04 2.32921745e-03]]\n"
     ]
    }
   ],
   "source": [
    "def forward(O,pi,a,b,N,M,T):\n",
    "    fwd=np.zeros((T,N))\n",
    "\n",
    "    #初始化\n",
    "    fwd[0]=pi*b[:,O[0]]\n",
    "\n",
    "    #递推:\n",
    "    for t in range(T-1):\n",
    "        fwd[t+1]=np.dot(fwd[t],a)*b[:,O[t+1]]\n",
    "    \n",
    "    return fwd\n",
    "\n",
    "print(forward(*hmm))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "前面计算了到最后一个时刻的某个状态的所有路径的概率累加和，为了计算所有路径的累加和，我们需要把最好一个时刻的所有概率再加起来："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.0033974411807013404\n"
     ]
    }
   ],
   "source": [
    "def full_prob(fwd):\n",
    "    return fwd[-1].sum()\n",
    "\n",
    "print(full_prob(forward(*hmm)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Viterbi算法\n",
    "\n",
    "最简单算法只考虑某一个时刻最可能的状态，也就是根据发射概率找出最可能的状态。$s_i=argmax_sP(o_i|s)$。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0, 0, 2, 0, 2])"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.argmax(ob,axis=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "上面的序列不会考虑上下文。但是HMM认为有些状态之间的跳转概率更大，因此即使某个发射概率可能比另外一个稍大一些，但是如果从前面的状态跳转到它的概率很低的话，我们也应该选择跳转概率更大的那个。所以我们需要考虑整个序列，综合考虑跳转概率和发射概率。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0 1 2 1 2]\n"
     ]
    }
   ],
   "source": [
    "def viterbi(O,pi,a,b,N,M,T):\n",
    "    d=np.zeros((T,N))\n",
    "    ph=np.zeros((T,N),dtype=np.int)\n",
    "\n",
    "    #初始化\n",
    "    for i in range(N):\n",
    "        d[0,i]=pi[i]*b[i,O[0]]\n",
    "        ph[0,i]=0\n",
    "\n",
    "    #递推\n",
    "    for t in range(1,T):\n",
    "        for j in range(N):\n",
    "            m=np.zeros(N)\n",
    "            for i in range(N):\n",
    "                m[i]=d[t-1,i]*a[i,j]\n",
    "            ph[t,j]=m.argmax()\n",
    "            d[t,j]=m.max()*b[j,O[t]]\n",
    "\n",
    "    #结束时刻\n",
    "    m=np.zeros(N)\n",
    "    for i in range(N):\n",
    "        m[i]=d[T-1,i]\n",
    "    Pv=m.max()\n",
    "\n",
    "    #back-tracking\n",
    "    Q=np.zeros(T,dtype=np.int)\n",
    "    Q[T-1]=m.argmax()\n",
    "    for t in reversed(range(T-1)):\n",
    "        Q[t]=ph[t+1,Q[t+1]]\n",
    "        \n",
    "    return Q\n",
    "\n",
    "print(viterbi(*hmm))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "稍快一点的Numpy实现："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0 1 2 1 2]\n"
     ]
    }
   ],
   "source": [
    "def viterbi(O,pi,a,b,N,M,T):\n",
    "    d=np.zeros((T,N))\n",
    "    ph=np.zeros((T,N),dtype=np.int)\n",
    "\n",
    "    #初始化\n",
    "    d[0]=pi*b[:,O[0]]\n",
    "    ph[0]=0\n",
    "\n",
    "    #递推\n",
    "    for t in range(1,T):\n",
    "        m=d[t-1]*a.T        \n",
    "        ph[t]=m.argmax(axis=1)\n",
    "        d[t]=m[np.arange(N),ph[t]]*b[:,O[t]]    \n",
    "\n",
    "    #结束时刻\n",
    "    Q=np.zeros(T,dtype=np.int)\n",
    "    Q[T-1]=np.argmax(d[T-1])\n",
    "    Pv=d[T-1,Q[T-1]]\n",
    "\n",
    "    #back-tracking\n",
    "    for t in reversed(range(T-1)):\n",
    "        Q[t]=ph[t+1,Q[t+1]]\n",
    "\n",
    "    return Q\n",
    "\n",
    "print(viterbi(*hmm))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "如果我们把先验概率和转移概率设置成均匀分布，那么Viterbi得到的路径和前面最大似然的路径是相同的。\n",
    "\n",
    "\n",
    "### 前向后向算法\n",
    "\n",
    "前向算法已经实现过了，这里我们实现后向算法："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[0.01161635 0.00942741 0.01065058]\n",
      " [0.0218539  0.04098234 0.03181981]\n",
      " [0.11128824 0.0895978  0.10222045]\n",
      " [0.20510646 0.39447244 0.30445164]\n",
      " [1.         1.         1.        ]]\n"
     ]
    }
   ],
   "source": [
    "def backward(O,pi,a,b,N,M,T):\n",
    "    bk=np.zeros((T,N))\n",
    "    #初始化\n",
    "    for i in range(N):\n",
    "        bk[T-1,i]=1\n",
    "    #推导:\n",
    "    for t in reversed(range(T-1)):\n",
    "        for i in range(N):\n",
    "            s=0\n",
    "            for j in range(N):\n",
    "                s+=a[i,j]*b[j,O[t+1]]*bk[t+1,j]            \n",
    "            bk[t,i]=s\n",
    "    return bk\n",
    "\n",
    "print(backward(*hmm))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Numpy的实现："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[0.01161635 0.00942741 0.01065058]\n",
      " [0.0218539  0.04098234 0.03181981]\n",
      " [0.11128824 0.0895978  0.10222045]\n",
      " [0.20510646 0.39447244 0.30445164]\n",
      " [1.         1.         1.        ]]\n"
     ]
    }
   ],
   "source": [
    "def backward(O,pi,a,b,N,M,T):\n",
    "    bk=np.zeros((T,N))\n",
    "    #初始化:\n",
    "    bk[T-1]=1\n",
    "    #递推:\n",
    "    for t in reversed(range(T-1)):\n",
    "        bk[t]=np.dot(bk[t+1]*b[:,O[t+1]],a.T)\n",
    "    return bk\n",
    "\n",
    "print(backward(*hmm))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "根据定义，前向概率是$P(o_1,...,o_t,s_t=j|M)$，而后向概率是$P(o_{t+1},...,o_T|s_t=j;M)$。状态占用概率$\\gamma_t(j)=P(s_t=j|O)=\\frac{\\alpha_t(j)\\beta_t(j)}{P(O)}$，而P(O)就是full_prob函数得到的结果。在语音识别章节有详细的推导过程，请读者参考相关内容。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[0.36016307 0.43808971 0.20174722]\n",
      " [0.2577077  0.5126293  0.229663  ]\n",
      " [0.27567901 0.06160171 0.66271928]\n",
      " [0.10042141 0.66054168 0.23903691]\n",
      " [0.27932224 0.03509784 0.68557992]]\n"
     ]
    }
   ],
   "source": [
    "def gamma(fwd,bk,fp):\n",
    "    gm=np.zeros((T,N))\n",
    "    for t in range(T):\n",
    "        for i in range(N):\n",
    "            gm[t,i]=(fwd[t,i]*bk[t,i])/fp\n",
    "    return gm\n",
    "\n",
    "print(gamma(forward(*hmm),backward(*hmm),full_prob(forward(*hmm))))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Numpy的一行实现："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[0.36016307 0.43808971 0.20174722]\n",
      " [0.2577077  0.5126293  0.229663  ]\n",
      " [0.27567901 0.06160171 0.66271928]\n",
      " [0.10042141 0.66054168 0.23903691]\n",
      " [0.27932224 0.03509784 0.68557992]]\n"
     ]
    }
   ],
   "source": [
    "def gamma(fwd,bk,fp):\n",
    "    return (fwd*bk)/fp\n",
    "\n",
    "print(gamma(forward(*hmm),backward(*hmm),full_prob(forward(*hmm))))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "我们可以把这3个概率绘制出来进行比较。我们可以发现前向概率随着时间的增加指数衰减(后向算法也类似，只是方向相反)。把它们加到一起就可以找到每个时刻最可能的状态。白色表示概率较大，因此从图中我们可以知道1-5个时刻最优的状态是[1,1,2,1,2]，这回Viterbi算法得到的[0,1,2,1,2]是很接近的。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5, 1.0, 'Full')"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA2QAAADSCAYAAAAlkwM9AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAEqBJREFUeJzt3X+wpXV9H/D3h12WXdglCGFQYY0Wra3TmWpmtdP6s2lIUTQyo5OGTmhIM6W/bHRqmiYzmWg6zV9t09RMxo5JDNNodFx/RSoJsc26xkYT+aURKIYQ6aK2QECWRQR2/fSPc1aPW8veu/fc+z333Ndr5gznnPuc7/187y6ffd7P832eW90dAAAANt5powsAAADYqgQyAACAQQQyAACAQQQyAACAQQQyAACAQQQyAACAQQQyFkpVvaKq7hldB7D4quqaqvq3C1BHV9WzR9cBbD5V9daqetf0+TOn/WT76LrYWALZkquqL1bVo1V1ZObx9NF1AVvLCb3owar6aFXtHV0XwLzY5+JUCWRbw2u6e/fM48sr/eB6HqVxBAi2nNd09+4kT0vyf5L88uB6VkXPAlbglPe52LoEsi2qqn6wqm6tqq9W1cer6q/OfO2LVfWvq+pzSR6pqn9UVdfOfP1Pq2r/zOtDVfX86fP/NH19uKpurKqXzmz31qp6f1W9q6oOJ7mqqnZNlx09WFW3JXnhhvwAgGG6++tJ3p/keUlSVZdV1c3TvnGoqt46u31VvaSq/nDarw5V1VUnjllVe6rqQFW9raqeNd32tOnXfrWq7p3Z9jer6k3T5z9WVbdX1cNVdVdV/eOZ7V5RVfdM++H/TvIb0/f/VVV9paq+XFX/cO4/IGBpfKdLMab7Wd8/qiYWj0C2BVXVX07yniRvSnJ+kuuSXFtVO2Y2uyLJZUnOSfLfk7y0qk6bnnrfkeRvTsf6S0l2J/nc9HOfSfL8JOcm+a0k+6tq58y4r81kR+ycJO9O8pYkF08ffzfJj857vsBiqaozk/y9JJ+evvVIkn+QSV+4LMk/rarLp9t+T5LfyeRs2vmZ9JdbThjvvEz61P/o7p/o7j9PcjjJC6abvCzJkZkDTy9PcnD6/N4kr05ydpIfS/Ifq+p7Z4Z/aib97HuSXF1Vlyb5ySSXJHlOEjtVAKyJQLY1fHh6tPirVfXhTHaEPtrdH+vuJ5L8+yS7kvytmc+8rbsPdfej3X1Xkocz2RF6WZLrk3y5qv5KJjs2f9Dd30iS7n5Xd/9Fdx/t7v+Q5Iwkz50Z91Pd/eHu/kZ3P5rkh5L8Qnc/0N2HkrxtXX8SwEgfrqqvJnkok0Dz75Kkuz/e3X8y7Qufy+SA0cunn/n7Sf5bd7+nu5+Y9pfZQPb0TMLV/u7+2Zn3DyZ5eVU9dfr6/dPXz8okfH12+r0/2t1/1hMHk/xekpfOjPONJG/p7sdmetZvdPfnu/uRJG+dz48GWBIn7nPBSQlkW8Pl3X3O9HF5Jjswdx//4jRMHUpy4cxnDp0wxsEkr8gkkB1M8vFMdphmjzSnqn5yuvznoemO13cl+e4nGffpJ7x3d4BldXl3n5NkZ5I3JDlYVU+tqr8xXW54X1U9lOSf5Ft9Y2+SP3uSMS/L5IDSfz7h/dme9Yl8e8/65kGkqnplVX26qh6Y9qxX5dt71n3TJZbH6VnAkzlxnwtOSiDbmr6cyfKbJElVVSY7PV+a2aZP+MzxnZuXTp8fzAmBbHq92E9lcgT5KdMdr4eS1JOM+5Xp9z7uGacyIWDz6O5j3f3BJMeSvCST5c0fSbK3u78rk3B1vG8cymRJ8//Pryb53STXVdVZM+8fzKRfvWL6/JNJXpxv71lnJPlAJqsELpj2rOuiZwHz80iSM4+/qKptmSy/hm8SyLam9yW5rKr+TlWdnuTNSR5L8odP8pmDSf52kl3dfU+SP0hyaZLzktw83WZPkqNJ7kuyvap+LpOlQSer5Weq6ilVdVGSf3GKcwI2iZp4bZKnJLk9k97xQHd/vapelMkyxePeneT7q+qHqmp7VZ13/CZCM96Q5I5MroXdlSTd/adJHk3yI0kOdvfhTO7s+Lp866z+jkyWVd+X5GhVvTLJD5yk/PdlckOi502vhXvLqfwMgC3jC0l2Tm9edHqSn82k78A3CWRbUHffkclOyi8nuT/JazK5TevjT/KZLyQ5kkkQy3Tn5q5MLqI/Nt3s+kyOVH8hk2U8X8//u0TxRD8/3fbPM7l24zdPbVbAJnBtVR3J5IYbv5DkR7v71iT/LMm/qaqHk/xcJqEnSdLd/yuTZYRvTvJAJjf0+Ouzg3Z3J7k6yT1JfnvmRkIHk/zF9PrU468ryU3Tzz2c5Cem3+/BTILgR55sAt39O0l+KcnvJ7lz+l+A76i7H8qkx/1aJiuRHsmkV8E31eTfMQAAADaaM2QAAACDCGQAAACDCGQAAACDCGQAAACDCGQAAACDbF+PQXft2tV79uxZj6E3xOHDh0eXwNRjjz02ugS+5f7u3tS/zLKqNv1tZSe/x31z27Zt2+gS1uyss846+UYL7tnPfvboEubixhtv3PS9aceOHb1z586Tb7jATjtt8x/jv/jiJ/sd9JvHXXfdNbqENVuWP4uV9qd1CWR79uzJ6173uvUYekMcOHBgdAlzcfTo0dElrNndd989uoQ1W5ZfLXHs2LHN/4exBLZvX5e2vaHOPffc0SWs2Qtf+MLRJazZtddeO7qEuaiqTd+bdu7cmX379o0uY02W4SDF/v37R5cwF1dcccXoEtbsQx/60OgS5mKl/WnzH84AAADYpAQyAACAQQQyAACAQQQyAACAQQQyAACAQQQyAACAQQQyAACAQQQyAACAQQQyAACAQQQyAACAQQQyAACAQQQyAACAQQQyAACAQQQyAACAQQQyAACAQQQyAACAQQQyAACAQQQyAACAQQQyAACAQQQyAACAQQQyAACAQU4ayKpqb1UdqKrbqurWqnrjRhQGAACw7LavYJujSd7c3TdV1Z4kN1bVx7r7tnWuDQAAYKmd9AxZd3+lu2+aPn84ye1JLlzvwgAAAJbdqq4hq6pnJnlBkj9aj2IAAAC2kpUsWUySVNXuJB9I8qbuPvwdvn51kquTZPfu3XMrEGAtZnsTwKKY7U1nnHHG4GqAkVZ0hqyqTs8kjL27uz/4nbbp7nd0977u3rdr16551ghwymZ70+haAI6b7U07duwYXQ4w0EruslhJfj3J7d39i+tfEgAAwNawkjNkL05yZZLvq6pbpo9XrXNdAAAAS++k15B19yeT1AbUAgAAsKWs6i6LAAAAzI9ABgAAMIhABgAAMIhABgAAMIhABgAAMIhABgAAMIhABgAAMIhABgAAMIhABgAAMIhABgAAMIhABgAAMIhABgAAMIhABgAAMIhABgAAMIhABgAAMIhABgAAMIhABgAAMIhABgAAMIhABgAAMIhABgAAMIhABgAAMIhABgAAMIhABgAAMMj29Rj0Gc94Rt7+9revx9Ab4pJLLhldwlwcOHBgdAlrtmPHjtElrNmjjz46ugSWyLZt20aXsGbHjh0bXcKaXXDBBaNLYImce+65ufLKK0eXsSZ79+4dXcKanXfeeaNLmIv7779/dAlrVlWjS9hQzpABAAAMIpABAAAMIpABAAAMIpABAAAMIpABAAAMIpABAAAMIpABAAAMIpABAAAMIpABAAAMIpABAAAMIpABAAAMIpABAAAMIpABAAAMIpABAAAMIpABAAAMIpABAAAMIpABAAAMIpABAAAMIpABAAAMIpABAAAMIpABAAAMIpABAAAMctJAVlXvrKp7q+rzG1EQAADAVrGSM2TXJLl0nesAAADYck4ayLr7E0ke2IBaAAAAtpS5XUNWVVdX1Q1VdcN99903r2EB1mS2N42uBeC42d505MiR0eUAA80tkHX3O7p7X3fvO//88+c1LMCazPam0bUAHDfbm3bv3j26HGAgd1kEAAAYRCADAAAYZCW3vX9Pkk8leW5V3VNVP77+ZQEAACy/7SfboLuv2IhCAAAAthpLFgEAAAYRyAAAAAYRyAAAAAYRyAAAAAYRyAAAAAYRyAAAAAYRyAAAAAYRyAAAAAYRyAAAAAYRyAAAAAYRyAAAAAYRyAAAAAYRyAAAAAYRyAAAAAYRyAAAAAYRyAAAAAYRyAAAAAYRyAAAAAYRyAAAAAYRyAAAAAYRyAAAAAYRyAAAAAbZvh6DPvjgg9m/f/96DL0hjh49OrqEuTjzzDNHl7BmjzzyyOgS1qyqRpcwF909ugSSPPHEE6NLWLMdO3aMLmHNbr755tElsEQef/zx3H333aPLWJPrr79+dAlr9vrXv350CXNx8cUXjy5hza666qrRJczFNddcs6LtnCEDAAAYRCADAAAYRCADAAAYRCADAAAYRCADAAAYRCADAAAYRCADAAAYRCADAAAYRCADAAAYRCADAAAYRCADAAAYRCADAAAYRCADAAAYRCADAAAYRCADAAAYRCADAAAYRCADAAAYRCADAAAYRCADAAAYRCADAAAYRCADAAAYRCADAAAYZEWBrKourao7qurOqvrp9S4KAABgKzhpIKuqbUl+JckrkzwvyRVV9bz1LgwAAGDZreQM2YuS3Nndd3X340nem+S161sWAADA8ltJILswyaGZ1/dM3wMAAGAN5nZTj6q6uqpuqKobDh8+PK9hAdZktjeNrgXguNne9LWvfW10OcBAKwlkX0qyd+b1RdP3vk13v6O793X3vrPPPnte9QGsyWxvGl0LwHGzvenMM88cXQ4w0EoC2WeSPKeqnlVVO5L8cJKPrG9ZAAAAy2/7yTbo7qNV9YYk1yfZluSd3X3rulcGAACw5E4ayJKku69Lct061wIAALClzO2mHgAAAKyOQAYAADCIQAYAADCIQAYAADCIQAYAADCIQAYAADCIQAYAADCIQAYAADCIQAYAADCIQAYAADCIQAYAADCIQAYAADCIQAYAADCIQAYAADCIQAYAADCIQAYAADCIQAYAADCIQAYAADCIQAYAADCIQAYAADCIQAYAADCIQAYAADBIdff8B616OMkdcx94Y313kvtHF7FG5rA4lmEez+3uPaOLWAu9aWEswxyS5ZjHMsxBb1ocy/D3yRwWwzLMIVlhf9q+Tt/8ju7et05jb4iqusEcxluGOSTLMY+qumF0DXOgNy2AZZhDshzzWJY5jK5hDjZ9b0qW5++TOYy3DHNIVt6fLFkEAAAYRCADAAAYZL0C2TvWadyNZA6LYRnmkCzHPMxhMZjD4liGeZjDYliGOSTLMQ9zWAzLMIdkhfNYl5t6AAAAcHKWLAIAAAwy10BWVZdW1R1VdWdV/fQ8x94oVfXOqrq3qj4/upZTVVV7q+pAVd1WVbdW1RtH17RaVbWzqv64qj47ncPPj67pVFXVtqq6uar+6+haTkVVfbGq/qSqbtmsdzPTmxaD3rRYNntvSvSnRbHZ+9My9KZEf1okq+1Nc1uyWFXbknwhySVJ7knymSRXdPdtc/kGG6SqXpbkSJL/0t1/bXQ9p6Kqnpbkad19U1XtSXJjkss3059FVVWSs7r7SFWdnuSTSd7Y3Z8eXNqqVdW/TLIvydnd/erR9axWVX0xyb7u3pS/D0RvWhx602LZ7L0p0Z8WxWbvT8vQmxL9aZGstjfN8wzZi5Lc2d13dffjSd6b5LVzHH9DdPcnkjwwuo616O6vdPdN0+cPJ7k9yYVjq1qdnjgyfXn69LHpLnisqouSXJbk10bXsoXpTQtCb1ocetPC0J8WwDL0pkR/2szmGcguTHJo5vU92YR/mZdNVT0zyQuS/NHYSlZverr6liT3JvlYd2+6OST5pSQ/leQbowtZg07ye1V1Y1VdPbqYU6A3LSC9abhl6E2J/sScbebelOhPC2RVvclNPZZYVe1O8oEkb+ruw6PrWa3uPtbdz09yUZIXVdWmWgZRVa9Ocm933zi6ljV6SXd/b5JXJvnn06UpcMr0prGWqDcl+hNztNl7U6I/LZBV9aZ5BrIvJdk78/qi6XsMMF07/IEk7+7uD46uZy26+6tJDiS5dHQtq/TiJD84XUf83iTfV1XvGlvS6nX3l6b/vTfJhzJZYrOZ6E0LRG9aCEvRmxL9iflZpt6U6E+jrbY3zTOQfSbJc6rqWVW1I8kPJ/nIHMdnhaYXdf56ktu7+xdH13Mqqur8qjpn+nxXJhc8/8+xVa1Od/9Md1/U3c/M5P+H3+/uHxlc1qpU1VnTC5xTVWcl+YEkm+0uWnrTgtCbFsMy9KZEf2J+lqE3JfrTojiV3jS3QNbdR5O8Icn1mVwM+b7uvnVe42+UqnpPkk8leW5V3VNVPz66plPw4iRXZnJU4Zbp41Wji1qlpyU5UFWfy+QfrI9196a89ekmd0GST1bVZ5P8cZKPdvfvDq5pVfSmhaI3MU/604JYgv60DL0p0Z8Wxap709xuew8AAMDquKkHAADAIAIZAADAIAIZAADAIAIZAADAIAIZAADAIAIZAADAIAIZAADAIAIZAADAIP8X484qLG9fUtwAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 1080x216 with 3 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig,ax=P.subplots(1,3,figsize=(15,3),sharey=True)\n",
    "\n",
    "\n",
    "ax[0].pcolormesh(forward(*hmm).T,cmap=P.cm.gray)\n",
    "ax[0].set_yticks(np.arange(0,N)+.5)\n",
    "ax[0].set_yticklabels(np.arange(0,N))\n",
    "ax[0].set_title('Forward')\n",
    "\n",
    "ax[1].pcolormesh(backward(*hmm).T,cmap=P.cm.gray)\n",
    "ax[1].set_title('Backward')\n",
    "\n",
    "g=gamma(forward(*hmm),backward(*hmm),full_prob(forward(*hmm)))\n",
    "ax[2].pcolormesh(g.T,cmap=P.cm.gray)\n",
    "ax[2].set_title('Full')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "注意我们可以抽取很多序列，它们的含义是不同的。\n",
    "\n",
    "  * 我们可以根据观察逐个时刻(独立，不考虑上下文)的计算每个时刻最可能的状态，这个状态下发射这个观察的概率最大，这不会考虑状态的跳转。我们把这个序列叫做最近观察(得到的状态)序列\n",
    "  * 我们也可以考虑最佳的前向或者后向的概率，它只考虑单个方向的信息。 \n",
    "  * 我们可以计算状态占用概率，得到某个时刻最可能的状态，注意这不是最优路径，它只是说明某个时刻的最可能的状态，但是它们串起来可能就不是最优的。 \n",
    "  * 最后我们可以通过Viterbi算法得到最优状态序列。\n",
    "  \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Best observation sequence: [0 0 2 0 2]\n",
      "Best forward prob. sequence: [1 1 2 1 2]\n",
      "Best backward prob. sequence: [0 1 0 1 0]\n",
      "Best full prob. sequence: [1 1 2 1 2]\n",
      "Best Viterbi sequence: [0 1 2 1 2]\n"
     ]
    }
   ],
   "source": [
    "print('Best observation sequence: {}'.format(np.argmax(ob,axis=1)))\n",
    "print('Best forward prob. sequence: {}'.format(np.argmax(forward(*hmm),axis=1)))\n",
    "print('Best backward prob. sequence: {}'.format(np.argmax(backward(*hmm),axis=1)))\n",
    "g=gamma(forward(*hmm),backward(*hmm),full_prob(forward(*hmm)))\n",
    "print('Best full prob. sequence: {}'.format(np.argmax(g,axis=1)))\n",
    "print('Best Viterbi sequence: {}'.format(viterbi(*hmm)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Baum-Welch算法\n",
    "\n",
    "现在我们可以确定序列的概率，这有助于帮助我们fine-tuning状态转移概率。读者可以阅读Rabiner论文里的Baum-Welch算法，但是这在现在主流的HMM-DNN混合系统中很少使用。\n",
    "\n",
    "总共有5个时刻，但是最后一个时刻没有跳转，而每个时刻都有3x3种可能的跳转，因此我们可以得到一个4x3x3的tensor。这就是我们介绍过$\\epsilon_t(i,j)$(这里用xi表示)，它表示t时刻的状态是i而t+1时刻是状态j的概率，它可以这样计算：\n",
    "$\\epsilon_t(i,j)=\\frac{\\alpha_t(i)a_{ij}b_j(o_{t+1})\\beta_{t+1}(j)}{P(O|\\lambda)}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[[0.05759556 0.26043679 0.04213072]\n",
      "  [0.19053119 0.11756602 0.1299925 ]\n",
      "  [0.00958095 0.1346265  0.05753978]]\n",
      "\n",
      " [[0.07858843 0.03352748 0.14559178]\n",
      "  [0.18399368 0.01071143 0.31792419]\n",
      "  [0.0130969  0.01736279 0.19920331]]\n",
      "\n",
      " [[0.0431881  0.20028435 0.03220656]\n",
      "  [0.02645703 0.01674273 0.01840194]\n",
      "  [0.03077628 0.44351459 0.18842841]]\n",
      "\n",
      " [[0.03108645 0.01036278 0.05897218]\n",
      "  [0.23466249 0.01067459 0.41520459]\n",
      "  [0.0135733  0.01406046 0.21140315]]]\n"
     ]
    }
   ],
   "source": [
    "def xi(fwd,bk,fp,O,pi,a,b,N,M,T):\n",
    "    ret=np.zeros((T-1,N,N))\n",
    "    for t in range(T-1):\n",
    "        for i in range(N):\n",
    "            for j in range(N):\n",
    "                ret[t,i,j]=(fwd[t,i]*a[i,j]*b[j,O[t+1]]*bk[t+1,j])/fp                \n",
    "    return ret\n",
    "\n",
    "print(xi(forward(*hmm),backward(*hmm),full_prob(forward(*hmm)),*hmm))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Numpy版本："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[[0.05759556 0.26043679 0.04213072]\n",
      "  [0.19053119 0.11756602 0.1299925 ]\n",
      "  [0.00958095 0.1346265  0.05753978]]\n",
      "\n",
      " [[0.07858843 0.03352748 0.14559178]\n",
      "  [0.18399368 0.01071143 0.31792419]\n",
      "  [0.0130969  0.01736279 0.19920331]]\n",
      "\n",
      " [[0.0431881  0.20028435 0.03220656]\n",
      "  [0.02645703 0.01674273 0.01840194]\n",
      "  [0.03077628 0.44351459 0.18842841]]\n",
      "\n",
      " [[0.03108645 0.01036278 0.05897218]\n",
      "  [0.23466249 0.01067459 0.41520459]\n",
      "  [0.0135733  0.01406046 0.21140315]]]\n"
     ]
    }
   ],
   "source": [
    "def xi(fwd,bk,fp,O,pi,a,b,N,M,T):\n",
    "    return fwd[:-1].reshape((T-1,N,1))*a.reshape((1,N,N))*b[:,O[1:]].T.reshape((T-1,1,N))*bk[1:].reshape((T-1,1,N))/fp\n",
    "\n",
    "print(xi(forward(*hmm),backward(*hmm),full_prob(forward(*hmm)),*hmm))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "我们可以估计先验概率pi——也就是初始状态的概率，它就是0时刻的状态占用概率："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0.36016307 0.43808971 0.20174722]\n"
     ]
    }
   ],
   "source": [
    "def exp_pi(gamma):\n",
    "    return gamma[0]\n",
    "\n",
    "print(exp_pi(gamma(forward(*hmm),backward(*hmm),full_prob(forward(*hmm)))))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "接下来我们可以估计跳转概率，它的计算公式为：\n",
    "$\\hat{a}_{ij}=\\frac{\\sum_t\\epsilon_t(i,j)}{\\sum_j \\sum_t\\epsilon_t(i,j)}=\\frac{\\sum_t\\epsilon_t(i,j)}{\\sum_t\\gamma_t(i)}$ 。\n",
    "公式看起来比较复杂，其实意思很简单：分子是所有的状态i->j的概率累加和，而分母是所有状态为i的概率和。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[0.21173506 0.50767206 0.28059288]\n",
      " [0.37997411 0.09307088 0.52695501]\n",
      " [0.05027686 0.45723049 0.49249264]]\n"
     ]
    }
   ],
   "source": [
    "def exp_a(gamma,xi,N):\n",
    "    e_a=np.zeros((N,N))\n",
    "    for i in range(N):\n",
    "        for j in range(N):\n",
    "            e_a[i,j]=np.sum(xi[:,i,j])/np.sum(gamma[:-1,i])\n",
    "            \n",
    "    return e_a\n",
    "\n",
    "fw=forward(*hmm)\n",
    "bk=backward(*hmm)\n",
    "fp=full_prob(fw)\n",
    "g=gamma(fw,bk,fp)\n",
    "x=xi(fw,bk,fp,*hmm)\n",
    "\n",
    "ea=exp_a(g,x,N)\n",
    "\n",
    "print(exp_a(g,x,N))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Numpy版本："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[0.21173506 0.50767206 0.28059288]\n",
      " [0.37997411 0.09307088 0.52695501]\n",
      " [0.05027686 0.45723049 0.49249264]]\n"
     ]
    }
   ],
   "source": [
    "def exp_a(gamma,xi,N):            \n",
    "    return xi[:].sum(axis=0)/gamma[:-1].sum(axis=0).reshape(N,1)\n",
    "\n",
    "fw=forward(*hmm)\n",
    "bk=backward(*hmm)\n",
    "fp=full_prob(fw)\n",
    "g=gamma(fw,bk,fp)\n",
    "x=xi(fw,bk,fp,*hmm)\n",
    "\n",
    "print(exp_a(g,x,N))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "最后，我们可以估计发射概率，它的计算公式为：\n",
    "$\\hat{b}_{jk}=\\frac{\\sum_t[o_t=k]\\gamma_t(j)}{\\sum_t\\gamma_t(j)}$\n",
    "它的意思是：观察为k并且状态为j的概率和除以状态为j的概率。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[0.21650863 0.21936989 0.56412148 0.        ]\n",
      " [0.03606741 0.02054956 0.94338302 0.        ]\n",
      " [0.32828259 0.33960677 0.33211064 0.        ]]\n"
     ]
    }
   ],
   "source": [
    "def exp_b(gamma,O,N,M):\n",
    "    e_b=np.zeros((N,M))\n",
    "    for j in range(N):\n",
    "        for k in range(M):        \n",
    "            e_b[j,k]=np.sum(gamma[O==k,j])/np.sum(gamma[:,j])\n",
    "    return e_b\n",
    "\n",
    "fw=forward(*hmm)\n",
    "bk=backward(*hmm)\n",
    "fp=full_prob(fw)\n",
    "g=gamma(fw,bk,fp)\n",
    "\n",
    "print(exp_b(g,O,N,M))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 更新\n",
    "\n",
    "Baum-Welch算法是一种EM算法，首先我们随机初始化一组参数，然后就可以用前向和后向算法计算前向后向概率，然后根据它们计算状态占用概率$\\gamma_t(i)$和状态跳转概率$\\epsilon_t(i,j)$，接下来就可以用这两个概率重新估计pi,$a_{ij},b_{jk}$。每次迭代都可以增加训练数据上的似然概率，我们可以不断的迭代直到收敛(到局部最优解)。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Initial probability: 0.0033974411807013404\n",
      "Update #1 probability: 0.011811485858807922 -- mean error: 0.042663336376065515\n",
      "Update #2 probability: 0.014185384539023913 -- mean error: 0.001425943000793851\n",
      "Update #3 probability: 0.020919176220825817 -- mean error: 0.00328568136827396\n",
      "Update #4 probability: 0.05353066527613736 -- mean error: 0.009594967056081836\n",
      "Update #5 probability: 0.16337070329542774 -- mean error: 0.010107459650387708\n",
      "Update #6 probability: 0.2379400536037555 -- mean error: 0.001483043805476551\n",
      "Update #7 probability: 0.24971058117425143 -- mean error: 4.5221029194295106e-05\n",
      "Update #8 probability: 0.24999983233192202 -- mean error: 2.7953543736447684e-08\n",
      "Update #9 probability: 0.24999999999994377 -- mean error: 9.389540929545353e-15\n",
      "Update #10 probability: 0.25 -- mean error: 1.0557709064664283e-27\n",
      "Update #11 probability: 0.25 -- mean error: 5.0180801380220883e-54\n",
      "Update #12 probability: 0.25 -- mean error: 8.023836839987856e-106\n",
      "Update #13 probability: 0.25 -- mean error: 2.0514979230799537e-209\n",
      "Update #14 probability: 0.25 -- mean error: 0.0\n",
      "Update #15 probability: 0.25 -- mean error: 0.0\n"
     ]
    }
   ],
   "source": [
    "print('Initial probability: {}'.format(full_prob(forward(*hmm))))\n",
    "\n",
    "hmm_new=hmm\n",
    "for i in range(15):\n",
    "    fw=forward(*hmm_new)\n",
    "    bk=backward(*hmm_new)\n",
    "    fp=full_prob(fw)\n",
    "    g=gamma(fw,bk,fp)\n",
    "    x=xi(fw,bk,fp,*hmm_new)\n",
    "\n",
    "    pi_new=exp_pi(g)\n",
    "    a_new=exp_a(g,x,N)\n",
    "    b_new=exp_b(g,O,N,M)\n",
    "    \n",
    "    err=np.concatenate(((pi_new-hmm_new[1]).ravel(),(a_new-hmm_new[2]).ravel(),(b_new-hmm_new[3]).ravel()))    \n",
    "    \n",
    "    hmm_new=(O,pi_new,a_new,b_new,N,M,T)\n",
    "\n",
    "    print('Update #{} probability: {} -- mean error: {}'.format(i+1,full_prob(forward(*hmm_new)),np.mean(err**2)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "通常我们用HMM来建模一个语音单元比如因子，通常使用3-5个状态来建模一个因子，通过发音词典，我们直到一个单词对应的因子序列，我们可以把每个因子的HMM拼接起来(所有通常在语音识别系统中HMM有一个开始状态和结束状态，这样便于拼接HMM)得到一个词的HMM。类似的我们也可以把一个句子中所有词的HMM拼接起来得到一个很大的句子的HMM。我们实际标注的粒度一般是句子级别的，因此我们可以用数据来训练这个很大的HMM的参数，但是这个大的HMM是\"临时\"拼凑的HMM，真正的参数还是在因子级别。\n",
    "\n",
    "而对于HMM-DNN的混合系统，我们通常不适用EM算法。对于深度神经网络，我们通常使用区分性(discriminative)训练。这要求我们的标注是帧级别的(也就是标注每一帧属于那个因子)，这显然不可能让人来标注，实际通常使用Viterbi算法来进行Force-Alignment。什么叫Force Alignment呢？普通的Viterbi算法的是在所有的路径里搜索最有的状态序列，而对于训练数据，我们是有句子级别的标注的，我们可以利用这个信息。\n",
    "\n",
    "举个例子，句子\"go home\"，如果不做限制，它是在所有可能的路径里搜索，可能(虽然理论上可能性较低)解码成\"/b/ae/t/\"，而Force Alignment已经知道句子是\"go home\"了，因此已经知道HMM模型是/g əʊ h əʊ m/。\n",
    "假设观察序列是：\n",
    "$o_1,o_2,...,o_10$\n",
    "那么所有可能的对齐方式为：\n",
    "```\n",
    "g g əʊ əʊ əʊ h əʊ əʊ m\n",
    "g əʊ əʊ əʊ əʊ h əʊ əʊ m\n",
    "...\n",
    "```\n",
    "它只需要从这些对齐里选择最可能的路径就行了。而不可能是：\n",
    "```\n",
    "b b ae ae ae ae ae ae ae t\n",
    "```\n",
    "\n",
    "\n",
    "或者使用其它损失函数(比如CTC)，完全抛弃HMM模型。\n",
    "\n",
    "## 用于语音识别的HMM拓扑结构\n",
    "\n",
    "如果要实际使用HMM来进行语音识别，我们需要扩展HMM模型的定义。我们之前定义的HMM只能用于识别一个词。\n",
    "\n",
    "假设我们创建了一个声学模型(acoustic model)而不是简单的发射概率矩阵，声学模型可以根据声学观察(比较短时傅里叶变换滤波器组的输出)识别出因子序列。我们可以用神经网络来实现声学模型，我们可以认为神经网络的输出就是发射概率(但实际并不是，神经网络是区分性模型，它的输出是$P(s|o_t)$而不是$P(o_t|s)$)。把它替换掉原来的GMM，我们就可以把观察序列转换(识别)成因子序列。但是我们需要识别词，因此我们还需要识别词，因此我们还需要一个模型把因子训练变成词序列(实际更加复杂，声学模型输出的通常是上下文相关的因子序列，我们首先需要把它转换成上下文无关的因子序列，然后通过发音词典把它转换成词序列，最后通过语言模型选择最可能的词序列。\n",
    "\n",
    "注意：声学模型输出的不是一个因子序列，而是很多可能的上下文相关的因子序列，当然每个序列有不同的概率得分；类似的一个上下文无关的因子序列可以转换成很多可能的词序列，它们也有不同的概率，这样的组合会非常大，我们需要从这些可能的组合中寻找概率最大的那个。\n",
    "\n",
    "但是我们发现把观察序列变成上下文相关的因子序列的HMM H、把上下文相关的因子序列变成上下文无关的C、把上下文无关的因子序列变成词的L和把词序列变成词序列的语言模型L。都可以统一的用WFST来表示(某些复杂的文法或者N-gram需要通过近似的方法变成WFST)，而它们的层级关系可以通过WFST的composition来实现。这样我们可以把一个复杂的多层的搜索网络用WFST来表示，而且我们可以使用WFST的Determinization和Minimization算法来压缩和提高搜索速度。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# WFST\n",
    "\n",
    "WFST定义了一种关系，给定一个输入字符串，它可以输出多个(或者零个)字符串已经对应的weight。\n",
    "\n",
    "### 例子\n",
    "\n",
    " \n",
    "我们这里使用[PyFst](http://pyfst.github.io/)库，它是对OpenFst的封装。我们可以定义如下的WFST："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [
    {
     "ename": "TypeError",
     "evalue": "cannot use a string pattern on a bytes-like object",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mTypeError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[0;32m~/py3.6-env/lib/python3.6/site-packages/IPython/core/formatters.py\u001b[0m in \u001b[0;36m__call__\u001b[0;34m(self, obj)\u001b[0m\n\u001b[1;32m    343\u001b[0m             \u001b[0mmethod\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mget_real_method\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mobj\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mprint_method\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    344\u001b[0m             \u001b[0;32mif\u001b[0m \u001b[0mmethod\u001b[0m \u001b[0;32mis\u001b[0m \u001b[0;32mnot\u001b[0m \u001b[0;32mNone\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 345\u001b[0;31m                 \u001b[0;32mreturn\u001b[0m \u001b[0mmethod\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    346\u001b[0m             \u001b[0;32mreturn\u001b[0m \u001b[0;32mNone\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    347\u001b[0m         \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32mfst/_fst.pyx\u001b[0m in \u001b[0;36mfst._fst._Fst._repr_svg_\u001b[0;34m()\u001b[0m\n",
      "\u001b[0;32mfst/_fst.pyx\u001b[0m in \u001b[0;36mfst._fst.StdVectorFst.draw\u001b[0;34m()\u001b[0m\n",
      "\u001b[0;32m~/py3.6-env/lib/python3.6/re.py\u001b[0m in \u001b[0;36msub\u001b[0;34m(pattern, repl, string, count, flags)\u001b[0m\n\u001b[1;32m    189\u001b[0m     \u001b[0ma\u001b[0m \u001b[0mcallable\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mit\u001b[0m\u001b[0;31m'\u001b[0m\u001b[0ms\u001b[0m \u001b[0mpassed\u001b[0m \u001b[0mthe\u001b[0m \u001b[0mmatch\u001b[0m \u001b[0mobject\u001b[0m \u001b[0;32mand\u001b[0m \u001b[0mmust\u001b[0m \u001b[0;32mreturn\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    190\u001b[0m     a replacement string to be used.\"\"\"\n\u001b[0;32m--> 191\u001b[0;31m     \u001b[0;32mreturn\u001b[0m \u001b[0m_compile\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mpattern\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mflags\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msub\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mrepl\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mstring\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcount\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    192\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    193\u001b[0m \u001b[0;32mdef\u001b[0m \u001b[0msubn\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mpattern\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mrepl\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mstring\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcount\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mflags\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mTypeError\u001b[0m: cannot use a string pattern on a bytes-like object"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<StdVectorFst with 1 states>"
      ]
     },
     "execution_count": 37,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "f=fst.Transducer()\n",
    "f.add_arc(0,0,'a','b')\n",
    "f.add_arc(0,0,'b','c')\n",
    "f.add_arc(0,0,'c','a')\n",
    "f[0].final=True\n",
    "f"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "下面我们再来创建另一个WFST，它的输入和输出是完全相同的。\n",
    "\n",
    "输入和输出完全相同的WFST其实就是一个WFSA，我们有时候也称之为*acceptor*.\n",
    "\n",
    "我们可以使用那还是*linear_chain*来生成一个类似链条的acceptor："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [
    {
     "ename": "TypeError",
     "evalue": "cannot use a string pattern on a bytes-like object",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mTypeError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[0;32m~/py3.6-env/lib/python3.6/site-packages/IPython/core/formatters.py\u001b[0m in \u001b[0;36m__call__\u001b[0;34m(self, obj)\u001b[0m\n\u001b[1;32m    343\u001b[0m             \u001b[0mmethod\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mget_real_method\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mobj\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mprint_method\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    344\u001b[0m             \u001b[0;32mif\u001b[0m \u001b[0mmethod\u001b[0m \u001b[0;32mis\u001b[0m \u001b[0;32mnot\u001b[0m \u001b[0;32mNone\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 345\u001b[0;31m                 \u001b[0;32mreturn\u001b[0m \u001b[0mmethod\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    346\u001b[0m             \u001b[0;32mreturn\u001b[0m \u001b[0;32mNone\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    347\u001b[0m         \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32mfst/_fst.pyx\u001b[0m in \u001b[0;36mfst._fst._Fst._repr_svg_\u001b[0;34m()\u001b[0m\n",
      "\u001b[0;32mfst/_fst.pyx\u001b[0m in \u001b[0;36mfst._fst.StdVectorFst.draw\u001b[0;34m()\u001b[0m\n",
      "\u001b[0;32m~/py3.6-env/lib/python3.6/re.py\u001b[0m in \u001b[0;36msub\u001b[0;34m(pattern, repl, string, count, flags)\u001b[0m\n\u001b[1;32m    189\u001b[0m     \u001b[0ma\u001b[0m \u001b[0mcallable\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mit\u001b[0m\u001b[0;31m'\u001b[0m\u001b[0ms\u001b[0m \u001b[0mpassed\u001b[0m \u001b[0mthe\u001b[0m \u001b[0mmatch\u001b[0m \u001b[0mobject\u001b[0m \u001b[0;32mand\u001b[0m \u001b[0mmust\u001b[0m \u001b[0;32mreturn\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    190\u001b[0m     a replacement string to be used.\"\"\"\n\u001b[0;32m--> 191\u001b[0;31m     \u001b[0;32mreturn\u001b[0m \u001b[0m_compile\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mpattern\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mflags\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msub\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mrepl\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mstring\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcount\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    192\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    193\u001b[0m \u001b[0;32mdef\u001b[0m \u001b[0msubn\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mpattern\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mrepl\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mstring\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcount\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mflags\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mTypeError\u001b[0m: cannot use a string pattern on a bytes-like object"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<StdVectorFst with 5 states>"
      ]
     },
     "execution_count": 38,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "i=fst.linear_chain(['a','a','c','b'],syms=f.isyms)\n",
    "i"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Composition\n",
    "\n",
    "WFST最有用的操作就是composition。这个操作的输入是两个WFST，第一个WFST的输入字母表是A，输出字母表是B；第二个的输入字母表是B输出字母表是C，则Composition之后的WFST的输入字母表是A而输出字母表是C。\n",
    "\n",
    "\n",
    "我们可以使用重载的>>运算符实现Composition操作，比如i>>f就是把i和f进行Composition。\n",
    "我们发现i>>f的作用就是把i的符号做一个替换——a变成b，b变成c，c变成a。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [
    {
     "ename": "TypeError",
     "evalue": "cannot use a string pattern on a bytes-like object",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mTypeError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[0;32m~/py3.6-env/lib/python3.6/site-packages/IPython/core/formatters.py\u001b[0m in \u001b[0;36m__call__\u001b[0;34m(self, obj)\u001b[0m\n\u001b[1;32m    343\u001b[0m             \u001b[0mmethod\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mget_real_method\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mobj\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mprint_method\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    344\u001b[0m             \u001b[0;32mif\u001b[0m \u001b[0mmethod\u001b[0m \u001b[0;32mis\u001b[0m \u001b[0;32mnot\u001b[0m \u001b[0;32mNone\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 345\u001b[0;31m                 \u001b[0;32mreturn\u001b[0m \u001b[0mmethod\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    346\u001b[0m             \u001b[0;32mreturn\u001b[0m \u001b[0;32mNone\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    347\u001b[0m         \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32mfst/_fst.pyx\u001b[0m in \u001b[0;36mfst._fst._Fst._repr_svg_\u001b[0;34m()\u001b[0m\n",
      "\u001b[0;32mfst/_fst.pyx\u001b[0m in \u001b[0;36mfst._fst.StdVectorFst.draw\u001b[0;34m()\u001b[0m\n",
      "\u001b[0;32m~/py3.6-env/lib/python3.6/re.py\u001b[0m in \u001b[0;36msub\u001b[0;34m(pattern, repl, string, count, flags)\u001b[0m\n\u001b[1;32m    189\u001b[0m     \u001b[0ma\u001b[0m \u001b[0mcallable\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mit\u001b[0m\u001b[0;31m'\u001b[0m\u001b[0ms\u001b[0m \u001b[0mpassed\u001b[0m \u001b[0mthe\u001b[0m \u001b[0mmatch\u001b[0m \u001b[0mobject\u001b[0m \u001b[0;32mand\u001b[0m \u001b[0mmust\u001b[0m \u001b[0;32mreturn\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    190\u001b[0m     a replacement string to be used.\"\"\"\n\u001b[0;32m--> 191\u001b[0;31m     \u001b[0;32mreturn\u001b[0m \u001b[0m_compile\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mpattern\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mflags\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msub\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mrepl\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mstring\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcount\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    192\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    193\u001b[0m \u001b[0;32mdef\u001b[0m \u001b[0msubn\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mpattern\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mrepl\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mstring\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcount\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mflags\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mTypeError\u001b[0m: cannot use a string pattern on a bytes-like object"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<StdVectorFst with 5 states>"
      ]
     },
     "execution_count": 39,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "o=i>>f\n",
    "o.project_output()\n",
    "o"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "我们也可以定义3个(链式的)WFST，然后使用union把它们并联成一个大的WFST，同时我们使用remove_epsilon去除$\\epsilon$。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [
    {
     "ename": "TypeError",
     "evalue": "cannot use a string pattern on a bytes-like object",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mTypeError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[0;32m~/py3.6-env/lib/python3.6/site-packages/IPython/core/formatters.py\u001b[0m in \u001b[0;36m__call__\u001b[0;34m(self, obj)\u001b[0m\n\u001b[1;32m    343\u001b[0m             \u001b[0mmethod\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mget_real_method\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mobj\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mprint_method\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    344\u001b[0m             \u001b[0;32mif\u001b[0m \u001b[0mmethod\u001b[0m \u001b[0;32mis\u001b[0m \u001b[0;32mnot\u001b[0m \u001b[0;32mNone\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 345\u001b[0;31m                 \u001b[0;32mreturn\u001b[0m \u001b[0mmethod\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    346\u001b[0m             \u001b[0;32mreturn\u001b[0m \u001b[0;32mNone\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    347\u001b[0m         \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32mfst/_fst.pyx\u001b[0m in \u001b[0;36mfst._fst._Fst._repr_svg_\u001b[0;34m()\u001b[0m\n",
      "\u001b[0;32mfst/_fst.pyx\u001b[0m in \u001b[0;36mfst._fst.StdVectorFst.draw\u001b[0;34m()\u001b[0m\n",
      "\u001b[0;32m~/py3.6-env/lib/python3.6/re.py\u001b[0m in \u001b[0;36msub\u001b[0;34m(pattern, repl, string, count, flags)\u001b[0m\n\u001b[1;32m    189\u001b[0m     \u001b[0ma\u001b[0m \u001b[0mcallable\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mit\u001b[0m\u001b[0;31m'\u001b[0m\u001b[0ms\u001b[0m \u001b[0mpassed\u001b[0m \u001b[0mthe\u001b[0m \u001b[0mmatch\u001b[0m \u001b[0mobject\u001b[0m \u001b[0;32mand\u001b[0m \u001b[0mmust\u001b[0m \u001b[0;32mreturn\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    190\u001b[0m     a replacement string to be used.\"\"\"\n\u001b[0;32m--> 191\u001b[0;31m     \u001b[0;32mreturn\u001b[0m \u001b[0m_compile\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mpattern\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mflags\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msub\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mrepl\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mstring\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcount\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    192\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    193\u001b[0m \u001b[0;32mdef\u001b[0m \u001b[0msubn\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mpattern\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mrepl\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mstring\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcount\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mflags\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mTypeError\u001b[0m: cannot use a string pattern on a bytes-like object"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<StdVectorFst with 16 states>"
      ]
     },
     "execution_count": 40,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "f=fst.Transducer()\n",
    "f.add_arc(0,0,'a','b')\n",
    "f.add_arc(0,0,'b','c')\n",
    "f.add_arc(0,0,'c','a')\n",
    "f.add_arc(0,0,'d','f')\n",
    "f.add_arc(0,0,'e','g')\n",
    "f.add_arc(0,0,'f','h')\n",
    "f.add_arc(0,0,'g','e')\n",
    "f.add_arc(0,0,'h','f')\n",
    "f.add_arc(0,0,'x','x')\n",
    "f.add_arc(0,0,'y','y')\n",
    "f.add_arc(0,0,'z','z')\n",
    "f[0].final=True\n",
    "\n",
    "i1=fst.linear_chain(['a','b','c','d','e'],syms=f.isyms)\n",
    "i2=fst.linear_chain(['a','b','x','y','z'],syms=f.isyms)\n",
    "i3=fst.linear_chain(['g','h','c','d','e'],syms=f.isyms)\n",
    "i=i1.union(i2)\n",
    "i=i.union(i3)\n",
    "i.remove_epsilon()\n",
    "i"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "我们可以发现里面有一些重复路径，比如从状态0->1和0->6都是a，因此这个WFST是非确定的WFST，它的搜索需要回溯。我们可以使用determinize函数先把它变成确定的WFST，然后使用minimize函数得到最小的WFST，从而压缩空间和提高搜索速度。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [
    {
     "ename": "TypeError",
     "evalue": "cannot use a string pattern on a bytes-like object",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mTypeError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[0;32m~/py3.6-env/lib/python3.6/site-packages/IPython/core/formatters.py\u001b[0m in \u001b[0;36m__call__\u001b[0;34m(self, obj)\u001b[0m\n\u001b[1;32m    343\u001b[0m             \u001b[0mmethod\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mget_real_method\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mobj\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mprint_method\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    344\u001b[0m             \u001b[0;32mif\u001b[0m \u001b[0mmethod\u001b[0m \u001b[0;32mis\u001b[0m \u001b[0;32mnot\u001b[0m \u001b[0;32mNone\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 345\u001b[0;31m                 \u001b[0;32mreturn\u001b[0m \u001b[0mmethod\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    346\u001b[0m             \u001b[0;32mreturn\u001b[0m \u001b[0;32mNone\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    347\u001b[0m         \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32mfst/_fst.pyx\u001b[0m in \u001b[0;36mfst._fst._Fst._repr_svg_\u001b[0;34m()\u001b[0m\n",
      "\u001b[0;32mfst/_fst.pyx\u001b[0m in \u001b[0;36mfst._fst.StdVectorFst.draw\u001b[0;34m()\u001b[0m\n",
      "\u001b[0;32m~/py3.6-env/lib/python3.6/re.py\u001b[0m in \u001b[0;36msub\u001b[0;34m(pattern, repl, string, count, flags)\u001b[0m\n\u001b[1;32m    189\u001b[0m     \u001b[0ma\u001b[0m \u001b[0mcallable\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mit\u001b[0m\u001b[0;31m'\u001b[0m\u001b[0ms\u001b[0m \u001b[0mpassed\u001b[0m \u001b[0mthe\u001b[0m \u001b[0mmatch\u001b[0m \u001b[0mobject\u001b[0m \u001b[0;32mand\u001b[0m \u001b[0mmust\u001b[0m \u001b[0;32mreturn\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    190\u001b[0m     a replacement string to be used.\"\"\"\n\u001b[0;32m--> 191\u001b[0;31m     \u001b[0;32mreturn\u001b[0m \u001b[0m_compile\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mpattern\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mflags\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msub\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mrepl\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mstring\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcount\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    192\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    193\u001b[0m \u001b[0;32mdef\u001b[0m \u001b[0msubn\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mpattern\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mrepl\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mstring\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcount\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mflags\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mTypeError\u001b[0m: cannot use a string pattern on a bytes-like object"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<StdVectorFst with 10 states>"
      ]
     },
     "execution_count": 41,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "i=i.determinize()\n",
    "i.minimize()\n",
    "i"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "对于WFST，我们可以通过project_output函数去掉输入，只保留输出，得到一个WFSA。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [
    {
     "ename": "TypeError",
     "evalue": "cannot use a string pattern on a bytes-like object",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mTypeError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[0;32m~/py3.6-env/lib/python3.6/site-packages/IPython/core/formatters.py\u001b[0m in \u001b[0;36m__call__\u001b[0;34m(self, obj)\u001b[0m\n\u001b[1;32m    343\u001b[0m             \u001b[0mmethod\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mget_real_method\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mobj\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mprint_method\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    344\u001b[0m             \u001b[0;32mif\u001b[0m \u001b[0mmethod\u001b[0m \u001b[0;32mis\u001b[0m \u001b[0;32mnot\u001b[0m \u001b[0;32mNone\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 345\u001b[0;31m                 \u001b[0;32mreturn\u001b[0m \u001b[0mmethod\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    346\u001b[0m             \u001b[0;32mreturn\u001b[0m \u001b[0;32mNone\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    347\u001b[0m         \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32mfst/_fst.pyx\u001b[0m in \u001b[0;36mfst._fst._Fst._repr_svg_\u001b[0;34m()\u001b[0m\n",
      "\u001b[0;32mfst/_fst.pyx\u001b[0m in \u001b[0;36mfst._fst.StdVectorFst.draw\u001b[0;34m()\u001b[0m\n",
      "\u001b[0;32m~/py3.6-env/lib/python3.6/re.py\u001b[0m in \u001b[0;36msub\u001b[0;34m(pattern, repl, string, count, flags)\u001b[0m\n\u001b[1;32m    189\u001b[0m     \u001b[0ma\u001b[0m \u001b[0mcallable\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mit\u001b[0m\u001b[0;31m'\u001b[0m\u001b[0ms\u001b[0m \u001b[0mpassed\u001b[0m \u001b[0mthe\u001b[0m \u001b[0mmatch\u001b[0m \u001b[0mobject\u001b[0m \u001b[0;32mand\u001b[0m \u001b[0mmust\u001b[0m \u001b[0;32mreturn\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    190\u001b[0m     a replacement string to be used.\"\"\"\n\u001b[0;32m--> 191\u001b[0;31m     \u001b[0;32mreturn\u001b[0m \u001b[0m_compile\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mpattern\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mflags\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msub\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mrepl\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mstring\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcount\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    192\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    193\u001b[0m \u001b[0;32mdef\u001b[0m \u001b[0msubn\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mpattern\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mrepl\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mstring\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcount\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mflags\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mTypeError\u001b[0m: cannot use a string pattern on a bytes-like object"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<StdVectorFst with 10 states>"
      ]
     },
     "execution_count": 42,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "o=i>>f\n",
    "o.project_output()\n",
    "o"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "我们可以用WFST来实现语音识别的解码器。比如我们可以使用一个WFST把观察序列转换成上下文相关的因子序列、接着使用一个WFST把它转换成上下文无关的因子序列、再用发音词典把它转化成词序列，最后用语言模型给这个词序列加上语言模型概率得分。我们把这几个WFST composite在一起，然后使用determinize和minimize优化这个WFST。有了这个WFSt之后，给定观察序列，我们直接就可以用WFST得到最可能的词序列。\n",
    "\n",
    "\n",
    "## HMM使用WFST的例子\n",
    "\n",
    "下面我们修改前面HMM的代码，首先我们构造把观察变成状态(因子)的WFSA："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [],
   "source": [
    "def obs_mat(O,b,N,T):\n",
    "    f=fst.StdAcceptor()\n",
    "        \n",
    "    for t in range(T):\n",
    "        for j in range(N):\n",
    "            f.add_arc(t,t+1,str(j),-np.log(b[j,O[t]]))\n",
    "        \n",
    "    f[T].final=True\n",
    "    \n",
    "    return f\n",
    "\n",
    "obs_fsa=obs_mat(O,b,N,T)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "注意：在WFST里weight的表示方法和HMM不同，这里我们使用[tropical semiring](https://en.wikipedia.org/wiki/Tropical_geometry)。也就是把(发射)概率用-log来表示。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [
    {
     "ename": "TypeError",
     "evalue": "cannot use a string pattern on a bytes-like object",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mTypeError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[0;32m~/py3.6-env/lib/python3.6/site-packages/IPython/core/formatters.py\u001b[0m in \u001b[0;36m__call__\u001b[0;34m(self, obj)\u001b[0m\n\u001b[1;32m    343\u001b[0m             \u001b[0mmethod\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mget_real_method\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mobj\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mprint_method\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    344\u001b[0m             \u001b[0;32mif\u001b[0m \u001b[0mmethod\u001b[0m \u001b[0;32mis\u001b[0m \u001b[0;32mnot\u001b[0m \u001b[0;32mNone\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 345\u001b[0;31m                 \u001b[0;32mreturn\u001b[0m \u001b[0mmethod\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    346\u001b[0m             \u001b[0;32mreturn\u001b[0m \u001b[0;32mNone\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    347\u001b[0m         \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32mfst/_fst.pyx\u001b[0m in \u001b[0;36mfst._fst._Fst._repr_svg_\u001b[0;34m()\u001b[0m\n",
      "\u001b[0;32mfst/_fst.pyx\u001b[0m in \u001b[0;36mfst._fst.StdVectorFst.draw\u001b[0;34m()\u001b[0m\n",
      "\u001b[0;32m~/py3.6-env/lib/python3.6/re.py\u001b[0m in \u001b[0;36msub\u001b[0;34m(pattern, repl, string, count, flags)\u001b[0m\n\u001b[1;32m    189\u001b[0m     \u001b[0ma\u001b[0m \u001b[0mcallable\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mit\u001b[0m\u001b[0;31m'\u001b[0m\u001b[0ms\u001b[0m \u001b[0mpassed\u001b[0m \u001b[0mthe\u001b[0m \u001b[0mmatch\u001b[0m \u001b[0mobject\u001b[0m \u001b[0;32mand\u001b[0m \u001b[0mmust\u001b[0m \u001b[0;32mreturn\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    190\u001b[0m     a replacement string to be used.\"\"\"\n\u001b[0;32m--> 191\u001b[0;31m     \u001b[0;32mreturn\u001b[0m \u001b[0m_compile\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mpattern\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mflags\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msub\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mrepl\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mstring\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcount\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    192\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    193\u001b[0m \u001b[0;32mdef\u001b[0m \u001b[0msubn\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mpattern\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mrepl\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mstring\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcount\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mflags\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mTypeError\u001b[0m: cannot use a string pattern on a bytes-like object"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<StdVectorFst with 6 states>"
      ]
     },
     "execution_count": 44,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "obs_fsa"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "这个FSA的输入是状态序列，输出这个状态序列的发射概率。比如给定当前的5个观察，我们可以计算状态(因子)序列为：\n",
    "```\n",
    "0 0 1 1 2\n",
    "```\n",
    "的发射概率为：\n",
    "```\n",
    "0.62 + 2.26 + 0.9 + 0.41 + 1.59\n",
    "```\n",
    "\n",
    "\n",
    "HMM除了发射概率之外还有状态的先验概率和跳转概率，我们也可以把它变成WFST："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [],
   "source": [
    "def trans_mat(pi,a,N):\n",
    "    f=fst.StdAcceptor()\n",
    "    \n",
    "    for i in range(N):\n",
    "        f.add_arc(0,i+1,str(i),-np.log(pi[i]))\n",
    "    \n",
    "    for i in range(N):\n",
    "        for j in range(N):\n",
    "            f.add_arc(i+1,j+1,str(j),-np.log(a[i,j]))\n",
    "        \n",
    "    for i in range(N):\n",
    "        f[i+1].final=True\n",
    "    \n",
    "    return f\n",
    "\n",
    "tr_fsa=trans_mat(pi,a,N)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [
    {
     "ename": "TypeError",
     "evalue": "cannot use a string pattern on a bytes-like object",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mTypeError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[0;32m~/py3.6-env/lib/python3.6/site-packages/IPython/core/formatters.py\u001b[0m in \u001b[0;36m__call__\u001b[0;34m(self, obj)\u001b[0m\n\u001b[1;32m    343\u001b[0m             \u001b[0mmethod\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mget_real_method\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mobj\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mprint_method\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    344\u001b[0m             \u001b[0;32mif\u001b[0m \u001b[0mmethod\u001b[0m \u001b[0;32mis\u001b[0m \u001b[0;32mnot\u001b[0m \u001b[0;32mNone\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 345\u001b[0;31m                 \u001b[0;32mreturn\u001b[0m \u001b[0mmethod\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    346\u001b[0m             \u001b[0;32mreturn\u001b[0m \u001b[0;32mNone\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    347\u001b[0m         \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32mfst/_fst.pyx\u001b[0m in \u001b[0;36mfst._fst._Fst._repr_svg_\u001b[0;34m()\u001b[0m\n",
      "\u001b[0;32mfst/_fst.pyx\u001b[0m in \u001b[0;36mfst._fst.StdVectorFst.draw\u001b[0;34m()\u001b[0m\n",
      "\u001b[0;32m~/py3.6-env/lib/python3.6/re.py\u001b[0m in \u001b[0;36msub\u001b[0;34m(pattern, repl, string, count, flags)\u001b[0m\n\u001b[1;32m    189\u001b[0m     \u001b[0ma\u001b[0m \u001b[0mcallable\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mit\u001b[0m\u001b[0;31m'\u001b[0m\u001b[0ms\u001b[0m \u001b[0mpassed\u001b[0m \u001b[0mthe\u001b[0m \u001b[0mmatch\u001b[0m \u001b[0mobject\u001b[0m \u001b[0;32mand\u001b[0m \u001b[0mmust\u001b[0m \u001b[0;32mreturn\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    190\u001b[0m     a replacement string to be used.\"\"\"\n\u001b[0;32m--> 191\u001b[0;31m     \u001b[0;32mreturn\u001b[0m \u001b[0m_compile\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mpattern\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mflags\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msub\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mrepl\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mstring\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcount\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    192\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    193\u001b[0m \u001b[0;32mdef\u001b[0m \u001b[0msubn\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mpattern\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mrepl\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mstring\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcount\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mflags\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mTypeError\u001b[0m: cannot use a string pattern on a bytes-like object"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<StdVectorFst with 4 states>"
      ]
     },
     "execution_count": 46,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "tr_fsa"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "我们可以把这两个FSA composition起来："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {},
   "outputs": [
    {
     "ename": "TypeError",
     "evalue": "cannot use a string pattern on a bytes-like object",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mTypeError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[0;32m~/py3.6-env/lib/python3.6/site-packages/IPython/core/formatters.py\u001b[0m in \u001b[0;36m__call__\u001b[0;34m(self, obj)\u001b[0m\n\u001b[1;32m    343\u001b[0m             \u001b[0mmethod\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mget_real_method\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mobj\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mprint_method\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    344\u001b[0m             \u001b[0;32mif\u001b[0m \u001b[0mmethod\u001b[0m \u001b[0;32mis\u001b[0m \u001b[0;32mnot\u001b[0m \u001b[0;32mNone\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 345\u001b[0;31m                 \u001b[0;32mreturn\u001b[0m \u001b[0mmethod\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    346\u001b[0m             \u001b[0;32mreturn\u001b[0m \u001b[0;32mNone\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    347\u001b[0m         \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32mfst/_fst.pyx\u001b[0m in \u001b[0;36mfst._fst._Fst._repr_svg_\u001b[0;34m()\u001b[0m\n",
      "\u001b[0;32mfst/_fst.pyx\u001b[0m in \u001b[0;36mfst._fst.StdVectorFst.draw\u001b[0;34m()\u001b[0m\n",
      "\u001b[0;32m~/py3.6-env/lib/python3.6/re.py\u001b[0m in \u001b[0;36msub\u001b[0;34m(pattern, repl, string, count, flags)\u001b[0m\n\u001b[1;32m    189\u001b[0m     \u001b[0ma\u001b[0m \u001b[0mcallable\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mit\u001b[0m\u001b[0;31m'\u001b[0m\u001b[0ms\u001b[0m \u001b[0mpassed\u001b[0m \u001b[0mthe\u001b[0m \u001b[0mmatch\u001b[0m \u001b[0mobject\u001b[0m \u001b[0;32mand\u001b[0m \u001b[0mmust\u001b[0m \u001b[0;32mreturn\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    190\u001b[0m     a replacement string to be used.\"\"\"\n\u001b[0;32m--> 191\u001b[0;31m     \u001b[0;32mreturn\u001b[0m \u001b[0m_compile\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mpattern\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mflags\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msub\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mrepl\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mstring\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcount\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    192\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    193\u001b[0m \u001b[0;32mdef\u001b[0m \u001b[0msubn\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mpattern\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mrepl\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mstring\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcount\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mflags\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mTypeError\u001b[0m: cannot use a string pattern on a bytes-like object"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<StdVectorFst with 16 states>"
      ]
     },
     "execution_count": 47,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "out_fsa=obs_fsa>>tr_fsa\n",
    "out_fsa"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "上面的WFST有些复杂，它的输入是状态序列，而输出是这个状态序列的概率(先验+跳转+发射，而且是-log空间上的)，因此我们可以使用*shortest_path*来寻找最优的路径(-log最小，也就是log最大，也就是概率最大的路径），这等价于Viterbi算法："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {},
   "outputs": [
    {
     "ename": "TypeError",
     "evalue": "cannot use a string pattern on a bytes-like object",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mTypeError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[0;32m~/py3.6-env/lib/python3.6/site-packages/IPython/core/formatters.py\u001b[0m in \u001b[0;36m__call__\u001b[0;34m(self, obj)\u001b[0m\n\u001b[1;32m    343\u001b[0m             \u001b[0mmethod\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mget_real_method\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mobj\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mprint_method\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    344\u001b[0m             \u001b[0;32mif\u001b[0m \u001b[0mmethod\u001b[0m \u001b[0;32mis\u001b[0m \u001b[0;32mnot\u001b[0m \u001b[0;32mNone\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 345\u001b[0;31m                 \u001b[0;32mreturn\u001b[0m \u001b[0mmethod\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    346\u001b[0m             \u001b[0;32mreturn\u001b[0m \u001b[0;32mNone\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    347\u001b[0m         \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32mfst/_fst.pyx\u001b[0m in \u001b[0;36mfst._fst._Fst._repr_svg_\u001b[0;34m()\u001b[0m\n",
      "\u001b[0;32mfst/_fst.pyx\u001b[0m in \u001b[0;36mfst._fst.StdVectorFst.draw\u001b[0;34m()\u001b[0m\n",
      "\u001b[0;32m~/py3.6-env/lib/python3.6/re.py\u001b[0m in \u001b[0;36msub\u001b[0;34m(pattern, repl, string, count, flags)\u001b[0m\n\u001b[1;32m    189\u001b[0m     \u001b[0ma\u001b[0m \u001b[0mcallable\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mit\u001b[0m\u001b[0;31m'\u001b[0m\u001b[0ms\u001b[0m \u001b[0mpassed\u001b[0m \u001b[0mthe\u001b[0m \u001b[0mmatch\u001b[0m \u001b[0mobject\u001b[0m \u001b[0;32mand\u001b[0m \u001b[0mmust\u001b[0m \u001b[0;32mreturn\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    190\u001b[0m     a replacement string to be used.\"\"\"\n\u001b[0;32m--> 191\u001b[0;31m     \u001b[0;32mreturn\u001b[0m \u001b[0m_compile\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mpattern\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mflags\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msub\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mrepl\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mstring\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcount\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    192\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    193\u001b[0m \u001b[0;32mdef\u001b[0m \u001b[0msubn\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mpattern\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mrepl\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mstring\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcount\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mflags\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mTypeError\u001b[0m: cannot use a string pattern on a bytes-like object"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<StdVectorFst with 6 states>"
      ]
     },
     "execution_count": 48,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sh_p=out_fsa.shortest_path()\n",
    "sh_p"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['0', '1', '2', '1', '2']\n"
     ]
    }
   ],
   "source": [
    "shp_l=[]\n",
    "for p in sh_p.paths():\n",
    "    for s in p:\n",
    "        shp_l.append(sh_p.osyms.find(s.olabel))\n",
    "print(shp_l)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0 1 2 1 2]\n"
     ]
    }
   ],
   "source": [
    "print(viterbi(*hmm))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "上面的演示了用WFST实现简单HMM解码的过程，实际的语音识别解码器需要H、C、L和G这4个WFST，实现从观察到词序列的解码过程。\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "py3.6-env",
   "language": "python",
   "name": "py3.6-env"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
